Non conservative Abelian sandpile model with BTW toppling rule 
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I. INTRODUCTION 



Recently Tsuchiya an Katori Q have introduced a non 
conservative Abelian sandpile model with a toppling rule 
similar to that of the well known Bak, Tang and Wiesen- 
feld (BTW) sandpile model §. The model is defined in 
a square lattice of size L and an integer energy profile 
Zij is considered. Sites with energy below the thresh- 
old z c = 4aC are stable. If the energy at any given 
site exceeds this threshold the site transfer energy 
to its four nearest neighbors following the toppling rule: 
Zij — > Zij — z c , z t ±ij — > Zi±\j +C and 2y±i — * %±i+C C 
is an integer number and a is such that z c is also integer. 
The boundaries are assumed to be open and the system 
is perturbed by adding a unit of energy at a site selected 
at random and letting it evolve until an equilibrium con- 
figuration is reached. 

On each toppling event an amount of energy e = 
4£(a — 1) > is dissipated. For a = 1 the model is 
conservative, it is just the BTW model but with a differ- 
ent scale of toppling and energy addition. In the BTW 
model 0] the energy added to perturb the system is 1 
and on toppling an active site transfers an energy 1 to 
each neighbor, they are of the same order. In the model 
defined above the energy added is still 1 but the energy 
transferred on toppling is £. In the limit a = 1 and 
C = 1 the BTW model m is recovered while for a = 1 
and £ 3> 1 it is similar to the BTW model but with a 
uniform driving. 

In the BTW model (a = 1 and £ — 1) the avalanches 
can be decomposed in a sequence of sub avalanches called 
waves H with well defined finite size scaling properties. 
On the contrary, the distribution of the overall avalanche 
size s is better described using a multi-fractal analysis 
H . The break down of the finite size scaling has been 
recently shown to be a consequence of the existence of 
correlations in the sequence of waves || . 

For a > 1 the model is non conservative but still 
Abelian Q. In the thermodynamic limit L — > 00, ex- 
act calculations by Tsuchiya and Katori yield the mean 
avalanche size (including avalanches with size zero) (T) = 



er 1 jjj. Since e = 4£(a — 1) they concluded that in the 
limit a — > 1 (T) diverges. However, as it is shown in sec. 
U this conclusion is wrong because a cannot goes to zero 
in an arbitrary way, in order to satisfy the constraint that 
z c = 4a£ remains integer. Here it is demonstrated that 
e > 1 and, therefore, (T) < 1 for all possible values of ( 
and a. 

The main goal of this work is to investigate the scaling 
properties of this non conservative BTW like model in the 
limit a — > 1 . The main questions are related with the ex- 
istence or not of criticality in the conservative limit a — > 1 
and if in this limit one recovers the conventional BTW 
model a = 1. From the analysis of the energy scales in- 
volved in the model ( sec. |n|) and numerical simulations 
(sec. Ill) it is concluded that the model is critical when 
a — > 1 but it does not belong to the universality class of 
the BTW model. Its relation with other non conserva- 
tive models with BTW like toppling rule introduced in 
the literature |^-|| is also discussed (sec. |l|). 



II. SCALING LAWS 

In this section some scaling laws are derived based on 
the energy scales involved in the model. The main idea 
of this approach is that the balance between input and 
dissipation of energy determine the scaling of some mag- 
nitudes with the dissipation per toppling event, follow- 
ing the general guidelines introduced by Vespiganani et 
al For this purpose the avalanches are assumed to 
be instantaneous and the analysis is focused in the time 
scale of the driving field. On each step one adds 1 unit of 
energy and measure the toppling activity and the energy 
dissipated. On each toppling event an amount of energy 
e = 4C(a — 1) is locally dissipated while an amount 4£ is 
transferred to nearest neighbors. For boundary sites part 
of the energy is also dissipated thr oug h the boundary. 

Let G(r) be the Green function po[ , the average num- 
ber of toppling events at a distance r from the site where 
the energy was added. Close to r = the effect of local 
dissipation gives an small contribution and the main en- 
ergy scale is given by the transport of the energy from 
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active sites to their nearest neighbors. On the contrary, 
far from r = the effects of the local dissipation be- 
comes more important. How far will depend on certain 
correlation length £, such that for r <C £ transport is 
more important than local dissipation while for r 3> £ 
the opposite occurs. 

Thus, there are two characteristic lengths in this 
model: the system size L and the correlation length £. 
The analysis developed above is valid in the thermody- 
namic limit L ;§> £. In this case the dissipation through 
the boundary of the system is negligible in comparison 
with the energy dissipated on each toppling event. In 
such a situation the only way to reach an stationary state 
is to balance the input of energy from the driving field 
with the energy locally dissipated. Moreover, since £ is 
the only characteristic length G(r) is expected to satisfy 
the scaling law 



G{r)=r^- d Hr/0, 



(1) 



where rj is an scaling exponent and d is the spatial di- 
mension. 

The amount of energy 6Ed(r) locally dissipated inside 
an hyper-circle of radius r is 



SE d (r) cx e / dpp^G{p) oc /(r/£), (2) 
Jo 

where f(x) = f~ dyy v ~ 1 J-(y)a,nd the second proportion- 
ality is obtained using eq. (Q). On the other hand, the 
average energy transported through its boundary 5E t (r) 
is given by 

SE t (r) cx Cr'-^Cr) « Cp-*g(r/Z), (3) 
dr 

where g(x) — d[x ri ~ 2 J-'(x)]/dx and the second propor- 
tionality is obtained using eq. (|l]) . The correlation length 
£ can be defined as the radius r at which this two contri- 
bution becomes of the same order. With this definition 
and equating eqs. (0) and fl3) with r — £ it results that 



1 



(4) 



On the other hand, on each step 1 unit of energy 
is added and in average the amount e(T) is dissipated, 
where (T) cx J °° drr d ~ 1 G(r) is the mean avalanche size 
counting even the case when it has size 0. Equating this 
two contributions it results that 



, > 1 1 



Moreover, using eq. one obtains (0) 
T] = 0. 



(5) 



(0) 



general arguments an can be easily adapted to any sand- 
pile model with local dissipation. The same argument 
(energy balance) has been previously used by Vespi- 
ganani et al H to understand the scaling properties of 
other sandpile models with local dissipation. Here, a 
slightly different approach has been considered where the 
new parameter £, the ratio between the energy received 
from nearest active neighbors and from the driving field, 
has been taken into account. 

From eq. @ Tsuchiya and Katori concluded that when 
a — ► 1 (T) diverges. However, this conclusion is not valid 
if z c — 4£a! is restricted to be an integer number. To 
show this let us write a = 1 + e/4£ which follows from 
eq. (||). But A(a — 4£ + e is restricted to take integer 
values. With £ being an integer number the only way 
to satisfy this requirement is that e is also integer, i.e. 
e = 1,2,3,.... Then, since the smaller non negative in- 
teger is 1 it is concluded that e > 1 and, therefore, from 
eq. (|) (T) < 1, i.e. it is bounded. 

Nevertheless, the correlation length £ in eq. (Q) does 
not only depend on e but also on £. For fixed e it di- 
verges in the limit £ — > oo and the model is critical. The 
real control parameter is then e e ff = e/£, i.e. the energy 
dissipated per toppling event relative to the characteris- 
tic energy scale of transport (. Although this result is 
in complete agreement with the field theory approach of 
Vespignani et al M the fact that (T) does not diverges 
when e e ff — * — > oo) excludes this model from their 
analysis. 

Moreover, in previous sandpile models conservation 
implies the scaling law (s) ~ £ 2 , where (s) is the mean 
avalanche size excluding those with size jj). To inves- 
tigate the validity of such scaling relation for the present 
model let us take into account that (s) is related to (T) 
through the expression 



(*) = (T)/Pa, 



(J) 



Eq. (g) reproduces the exact result by Tsuchiya and 
Katori. The present approach is however based on more 



where P a is the probability to obtain an avalanche with 
non zero size. In the models considered by Vespignani et 
a l @ C = 1 an d, therefore, from eqs. (||), (||) and ([?]) it 
results that (s) ~ £ 2 /P Q . Moreover, in this models P a has 
a finite value and, therefore, one obtains the mentioned 
scaling law (s) ~ £ 2 . 

On the contrary, in the model considered here (s) can 
not be related to £ using these arguments. For fixed 
e from eqs. (||) and (0) one obtains that (s) ~ l/P a - 
Thus, from the energy balance invoked above we cannot 
say anything about the scaling of (s) with £ (an exponent 
2 will be an accidental coincidence) and, therefore, this 
model belongs to a new universality class. 



III. NUMERICAL SIMULATIONS AND 
DISCUSSION 

In this sections results obtained from numerical sim- 
ulations of the model studied above are presented. The 
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simulations were performed using e = 1, L = 4096 and 
£ ranging from 2° to 2 10 (e e // = 1/C ranging from 1 to 
2~ 10 ). For these values the condition L <C £ was observed 
to be satisfied. Statistics was taken over 10 8 avalanches 
after the system reached the stationary state. 

Before entering in the analysis of the statistics of the 
avalanches let us check the validity of eq. (||) . The log- log 
plot of (s) vs. £ is shown in fig. A clear linear behav- 
ior is observed for log 10 £ > 5 suggesting that above this 
value simple scaling applies. On top of this points the 
numerically computed values of 1 / P a are plotted obtain- 
ing an overlap in agreement with eq. (J|). If the scaling 
relation (s) ~ £ 2 where valid, using eq. (^), (s) ~ (. 
However, a linear fit to this log-log plot gives a slope 
- 0.9. 

The fact that this scaling relation thus not hold is 
clearly shown in fig. |^, where the stationary energy dis- 
tribution is shown. As it can be seen P a = P Zc .-i does 
not decreases as 1/C but slower, which explains why (s) 
increases with £ with a slope smaller than 1. The rest 
of the distribution scales like 1/C which is just a conse- 
quence of the increase of the density of possible values of 
z. 

The avalanche statistics will be characterized by: the 
number of toppling events s and steps t required to reach 
an stable configuration, the number of sites a "touched" 
by the avalanche, and the characteristic radius of the 
cluster formed by these sites r. The main goal of the sim- 
ulations is to determine the probability densities p x {x, C) 
[x = s, t, a, r) in the stationary state. 

One can easily see that s = a, in other words sites 
topples only once within an avalanche. In this model, 
as a difference with the original BTW model, only one 
wave of topplings takes place. The first wave if generated 
from an initial site with height z = z c = 4( + e. When 
this site topples it transfers an amount equal to z c to its 
nearest neighbors and, therefore, ends with energy z = 0. 
The best we can have to obtain a second toppling at this 
site is that its four nearest neighbors also become active. 
In such a case the initial side will receive 4£ < z c units 
of energy, which is not enough to make it active again. 
Hence, no second wave will be obtained yielding s = a. 

Since the waves are known to satisfy well defined 
finite-size scaling properties and in the present model 
an avalanche is made by one wave it is expected that 
the distributions p x (x, C) also satisfy a finite-size scaling. 
However, the scaling exponents will not necessarily be 
those obtained for the scaling of waves because, in the 
present model, conservation does not introduce any scal- 
ing relation among the different scaling exponents. 

If finite size scaling applies then these densities will 
satisfy 



Px{z,() = x T *g[x/x c (()], 



(8) 



decay given by Q. The validity of this scaling form is 
supported by the numerical results. The cut off x c is 
determined by the existence of the characteristic length 
£ ~ and is expected to scale as x c ~ £ Dx ~ £ da; , where 
d x = D x v is an effective fractal dimension. 

To compute the exponents t x and d x the moment anal- 
ysis technique introduced by De Menech et al jll]] is used. 
The moments of the probability density in eq. (||) are 
given by 



where 



M) 



dxp(x)x q ~ C : 



c)d x 



Xq) 



(9) 



(10) 



where t x is the power law exponent characterizing the 
self-similar regime and x c is a cut off above which the 
distribution deviates from a power law and has a fast 



The last equivalence in eq. (^|) is valid for values of q not 
to small, for which the precise form of p x (z,C) a ^ small 
x is not important. 

cr x (q) can be determined from a linear fit to the log- 
log plot of (x q ) vs. C- The resulting values using 
£ = 2 5 , 2 6 ,..., 2 10 are shown in fig. [| In all cases 
[x = s,t,r) for q larger than 1 a well defined linear de- 
pendence is observed. From the linear fit (see eq. 1C]) to 
these plots the exponents t x and d x are computed. The 
results are shown in table |[ 

The exponent v is very close to 1/2 in very good 
agreement with the scaling arguments of previous sec- 
tion. On the other hand, d s is quite close to 1 which 
implies that the avalanche size (or area) scale as s ~ r 2 , 
i.e. avalanches are compact D s — 2. With this value, the 
scaling relation (2 — t s )D s = 2 yields the power law ex- 
ponent t s = 1 which is clearly in disagreement with the 
value computed numerically. The reason for this result is 
that conservation does not introduce any scaling relation 
as it generally occurs in sandpile models Q . 

The exponents computed using the moment analysis 
technique can be checked using rescaled plots of the in- 
tegrated distribution P x (x, Q =( a °dxp x {x,C). The re- 
sulting plots are shown in figs. and The scaling 
works quite good supporting the validity of the reported 
exponents. 

In the literature we can find other sandpile models with 
local dissipation and BTW like toppling rule ■ In the 
models considered in and Q the energy profile is con- 
tinuous and the dissipation rate per toppling event e is 
a control parameter which can take any real value and, 
therefore, can be tuned to zero. Another feature of these 
models is that only one wave of toppling can take place 
and, therefore, for any finite e the model is in a different 
universality class from that of the BTW model. 

On the other hand in || the energy profile is discrete 
as in the original BTW model at the prize of introducing 
stochasticity in the model. In this case with a probability 
p energy is fully dissipated yielding an average dissipa- 
tion per toppling event e = 2dp. Clearly p may take 
any real variable between and 1 and, therefore, also in 
this case the dissipation per toppling event can be fine 
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tuned to zero. As a difference with the models described 
in the previous paragraph, in this case multiple toppling 
of a site within an avalanche is possible, which make it 
closer to the original BTW model. Moreover, the use of 
finite size scaling techniques can be also questioned and 
a multi-fractal analysis may be more appropriate [fj"3| , 
which is another characteristic feature of the BTW model 
j| . All this elements together with the numerical results 
reported in || suggest that in the limit p — > 1 (e — > 0) 
this model belong to the same universality class of the 
BTW model. 

A common feature of all this models ^-^] is that 
(s) ~ (T 1 as predicted by the field theory approach 
of Vespignani et al |J, leading to the scaling relation 
(2 — t s )D s = 2. On the contrary, in the present model 
the scaling of (s) with e e // is not known and conserva- 
tion does not introduce the above scaling relation. Hence, 
the model introduced by Tsuchiya and Katori belongs to 
different class among sandpile models. 
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IV. SUMMARY AND CONCLUSIONS 

A non conservative Abelian sandpile model with a 
BTW like toppling rule has been studied. The model 
can be though as the only possible generalization of the 
BTW model to include local dissipation without intro- 
ducing stochasticity in the toppling rule and keeping a 
discrete energy profile. However, the scaling approach 
and numerical the simulations reported here show that 
it does not belong to the universality class of the BTW 
model, not even to the universality class of any sandpile 
model previously considered in the literature. 
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FIG. 1. Log-log plot of the mean avalanche size (excluding 
avalanches with size s = 0) as a function £. It can be clearly 
seen that it scales as P^ 1 , the probability per unit step to 
obtain an avalanche with s > 0. The line is a linear fit to the 
high £ interval. 
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FIG. 2. Probability P z that a site has energy z in the sta- 
tionary state, for different values of £ = 2 n . z is expressed in 
units of the threshold z c = 4£+ 1 while P z has being rescaled 
by an amount £ because with increasing £ the density of z/z c 
values increases as 
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FIG. 3. Moment exponent cr x (q) for different values 
of q and x = s,t,r. The lines are linear fits 
[& x (q) — (1 — T x )d x + d x q] to the interval 1 > q < 3. The 
resulting exponents and d x are shown in tab. 
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FIG. 5. Scaled plot of the integrated distribution of 
avalanche durations using the exponents displayed in tab ffl. 
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FIG. 4. Scaled plot of the integrated distribution of 
avalanche sizes (or area since s = a in this model) using the 
exponents displayed in tab |. 
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FIG. 6. Scaled plot of the integrated distribution of 
avalanche radius using the exponents displayed in tab 



d s 


dt 


d r = V T s T t T r 


0.994(5) 


0.630(5) 


0.495(5) 1.11(1) 1.16(1) 1.14(1) 


D s = d s /u 


z = D t = d t /v 




2.01(1) 


1.27(1) 





TABLE I. Scaling exponents obtained from linear fits 
[o x (q) = (1 — r x )d x + d x q] to the data shown in fig. ^. 
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